home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / PROGRAMR / GSRC208A.ZIP / KINETICS.C < prev    next >
C/C++ Source or Header  |  1993-07-18  |  28KB  |  924 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*      reaction rate equations      */
  11. /*                                   */
  12. /*        Zortech C/C++ 3.0 r4       */
  13. /*          MICROSOFT C 6.00         */
  14. /*          Visual C/C++ 1.0         */
  15. /*           QuickC/WIN 1.0          */
  16. /*             ULTRIX cc             */
  17. /*              GNU gcc              */
  18. /*                                   */
  19. /*   (include here compilers that    */
  20. /*   compiled GEPASI successfully)   */
  21. /*                                   */
  22. /*************************************/
  23.  
  24. #include <math.h>
  25. #include "globals.h"
  26. #include "globvar.h"
  27. #include "datab.h"
  28.  
  29. /*****************************************/
  30. /*                                       */
  31. /*              CONVENTIONS              */
  32. /*                                       */
  33. /* all functions take as arguments:      */
  34. /* double s[] - vector of concentrations */
  35. /* int r - the index of the reaction     */
  36. /*                                       */
  37. /*     all fuctions return a double.     */
  38. /*                                       */
  39. /*                                       */
  40. /*****************************************/
  41.  
  42.  
  43. /*                                   */
  44. /* user defined rate-equations         */
  45. /*                                   */
  46.  
  47. double value_of( double s[], int r, struct treet *t, int idx )
  48. {
  49.  switch( t->node[idx].item )
  50.  {
  51.   case 'N': return (double) t->constant[(int)t->node[idx].val];
  52.   case 'I': switch( (int) t->id[(int)t->node[idx].val][9] )
  53.             {
  54.              case 0: return params[r][ (int) t->id[(int)t->node[idx].val][0] ];
  55.              case 1:
  56.              case 2:
  57.              case 3: return s[ eff[r][ (int) t->id[(int)t->node[idx].val][0] ] ];
  58.             }
  59.   case 'F': switch( t->node[idx].val )
  60.             {
  61.              case '+': return value_of( s, r, t, (int) t->node[idx].left );
  62.              case '-': return - value_of( s, r, t, (int) t->node[idx].left );
  63.              case 'L': return log10( value_of( s, r, t, (int) t->node[idx].left ) );
  64.              case 'l': return log( value_of( s, r, t, (int) t->node[idx].left ) );
  65.              case 'e': return exp( value_of( s, r, t, (int) t->node[idx].left ) );
  66.              case 'S': return sin( value_of( s, r, t, (int) t->node[idx].left ) );
  67.              case 'C': return cos( value_of( s, r, t, (int) t->node[idx].left ) );
  68.             }
  69.   case 'O': switch( t->node[idx].val )
  70.             {
  71.              case '+': return value_of( s, r, t, (int) t->node[idx].left ) +
  72.                               value_of( s, r, t, (int) t->node[idx].right );
  73.              case '-': return value_of( s, r, t, (int) t->node[idx].left ) -
  74.                               value_of( s, r, t, (int) t->node[idx].right );
  75.              case '*': return value_of( s, r, t, (int) t->node[idx].left ) *
  76.                               value_of( s, r, t, (int) t->node[idx].right );
  77.              case '/': return value_of( s, r, t, (int) t->node[idx].left ) /
  78.                               value_of( s, r, t, (int) t->node[idx].right );
  79.              case '^': return pow( value_of( s, r, t, (int) t->node[idx].left ),
  80.                                    value_of( s, r, t, (int) t->node[idx].right ) );
  81.             }
  82.  }
  83. }
  84.  
  85. double translate(double s[], int r)
  86. {
  87.  return value_of( s, r, &tree[kinetype[r]-MAX_TYP], (int) tree[kinetype[r]-MAX_TYP].node[0].left );
  88. }
  89.  
  90.  
  91. /*                                   */
  92. /* the zero returning function       */
  93. /*                                   */
  94.  
  95. double ret_zero( double s[], int r, int e )
  96. {
  97.  return (double) 0;
  98. }
  99.  
  100.  
  101. /*                                   */
  102. /* derivatives for effectors         */
  103. /* by finite differences             */
  104. /*                                   */
  105.  
  106. double deff( double s[], int r, int e )
  107. {
  108.  double k, x, x1, f1, f2;
  109.  
  110.  x = k = s[e];
  111.  /*
  112.     if x is zero, we will calculate the derivative at a small positive
  113.     value (no point in considering negative values!). let's stick with
  114.     hrcz (the flux resolution)
  115.  */
  116.  if( x==0 ) x = options.hrcz;
  117.  x1 = x * 0.01;
  118.  s[e] = x - x1;
  119.  f1 = rateq[r]( s, r );
  120.  s[e] = x + x1;
  121.  f2 = (rateq[r]( s, r ) - f1)/(x*0.02);
  122.  s[e] = k;
  123.  return f2;
  124. }
  125.  
  126. /*                                   */
  127. /* A -->  or  --> B                  */
  128. /*                                   */
  129.  
  130. double i0(double s[], int r)
  131. {
  132.  return params[r][0];
  133. }
  134.  
  135.  
  136. /*                                   */
  137. /* A --> B                           */
  138. /*                                   */
  139.  
  140. double i11(double s[], int r)
  141. {
  142.  return s[eff[r][0]]*params[r][0];
  143. }
  144.  
  145. /*                                 */
  146. /* A <--> B                        */
  147. /*                                 */
  148.  
  149. double r11(double s[], int r)
  150. {
  151.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*params[r][1];
  152. }
  153.  
  154. /*                                 */
  155. /* A + B --> C                     */
  156. /*                                 */
  157.  
  158. double i21(double s[], int r)
  159. {
  160.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  161. }
  162.  
  163. /*                                 */
  164. /* A + B <--> C                    */
  165. /*                                 */
  166.  
  167. double r21(double s[], int r)
  168. {
  169.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*params[r][1];
  170. }
  171.  
  172. /*                                   */
  173. /* A --> B + C                       */
  174. /*                                   */
  175.  
  176. double i12(double s[], int r)
  177. {
  178.  return s[eff[r][0]]*params[r][0];
  179. }
  180.  
  181. /*                                 */
  182. /* A <--> B + C                    */
  183. /*                                 */
  184.  
  185. double r12(double s[], int r)
  186. {
  187.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*params[r][1];
  188. }
  189.  
  190. /*                                 */
  191. /* A + B + C --> D                 */
  192. /*                                 */
  193.  
  194. double i31(double s[], int r)
  195. {
  196.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
  197. }
  198.  
  199. /*                                 */
  200. /* A + B + C <--> D                */
  201. /*                                 */
  202.  
  203. double r31(double s[], int r)
  204. {
  205.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*params[r][1];
  206. }
  207.  
  208. /*                                   */
  209. /* A --> B + C + D                   */
  210. /*                                   */
  211.  
  212. double i13(double s[], int r)
  213. {
  214.  return s[eff[r][0]]*params[r][0];
  215. }
  216.  
  217. /*                                 */
  218. /* A <--> B + C + D                */
  219. /*                                 */
  220.  
  221. double r13(double s[], int r)
  222. {
  223.  return s[eff[r][0]]*params[r][0]-s[eff[r][1]]*s[eff[r][2]]*s[eff[r][3]]*params[r][1];
  224. }
  225.  
  226. /*                                 */
  227. /* A + B --> C + D                 */
  228. /*                                 */
  229.  
  230. double i22(double s[], int r)
  231. {
  232.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  233. }
  234.  
  235. /*                                 */
  236. /* A + B <--> C + D                */
  237. /*                                 */
  238.  
  239. double r22(double s[], int r)
  240. {
  241.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0]-s[eff[r][2]]*s[eff[r][3]]*params[r][1];
  242. }
  243.  
  244. /*                                 */
  245. /* A + B + C --> D + E             */
  246. /*                                 */
  247.  
  248. double i32(double s[], int r)
  249. {
  250.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0];
  251. }
  252.  
  253. /*                                 */
  254. /* A + B + C <--> D + E            */
  255. /*                                 */
  256.  
  257. double r32(double s[], int r)
  258. {
  259.  return s[eff[r][0]]*s[eff[r][1]]*s[eff[r][2]]*params[r][0]-s[eff[r][3]]*s[eff[r][4]]*params[r][1];
  260. }
  261.  
  262. /*                                 */
  263. /* A + B --> C + D + E             */
  264. /*                                 */
  265.  
  266. double i23(double s[], int r)
  267. {
  268.  return s[eff[r][0]]*s[eff[r][1]]*params[r][0];
  269. }
  270.  
  271. /*                                 */
  272. /* A + B <--> C + D + E            */
  273. /*                                 */
  274.  
  275. double r23(double s[], int r)
  276. {
  277.  return s[eff[r][0]]*s[eff[r][1]]*